%%% Written by  Daniel Lewis, 2020, for "Robust Inference in Models 
%%% Identified via Heteroskedasticity", Review of Economics and Statistics.
%
% This file reads loads the empirical dataset, conducts F-stat based tests
% for weak identification, and estimates the four specifications reported
% in the paper. It computes both standard and non-robust confidence sets
% for each parameter, and thus also provides results for the Andrews (2018)
% test for weak identification. ?It compiles these results into Tables 4
% and 5 of the paper.

clear; clc

%% Load the data

load('dataset')
T=length(HF);
Tc=length(HFc);
Tp=T-Tc;

%% First-stage F-statistics

Table4=NaN(4,8);

% Daily yield changes
NYd=[NY(1:Tc,:)-mean(NY(1:Tc,:));NY(Tc+1:end,:)-mean(NY(Tc+1:end,:))]; % demean the data (the model assumes eta is mean zero in each sub-sample)
Z=[repmat(-T/Tc,Tc,1);repmat(T/Tp,Tp,1)].*NYd; % construct instrument
[b,~,r]=regress(NYd,Z); % first stage
varb=T*mean((Z.*r).^2)/sum(Z.^2)^2;  % heteroskedasticity-robust asymptotic variance
Fstat=b^2/varb; % F-statistic
fprintf('Daily changes F-stat = %.4f\n\n',Fstat);
Table4(1:2,1)=Fstat;

% NS (2018) 30-min "Policy News"
HFd=[HF(1:Tc,:)-mean(HF(1:Tc,:));HF(Tc+1:end,:)-mean(HF(Tc+1:end,:))]; % demean the data (the model assumes eta is mean zero in each sub-sample)
Z=[repmat(-T/Tc,Tc,1);repmat(T/Tp,Tp,1)].*HFd; % construct instrument
[b,~,r]=regress(HFd,Z); % first stage
varb=T*mean((Z.*r).^2)/sum(Z.^2)^2; % heteroskedasticity-robust asymptotic variance
Fstat=b^2/varb; % F-statistic
fprintf('30 min. changes F-stat = %.4f\n\n',Fstat);
Table4(3:4,1)=Fstat;

Table4(:,2:4)=Table4(:,1)>[15.06 23.11 37.42];

%% Set-up

% Compute Jacobian symbolically and construct function:
Hs=sym('H',2); % sympolic expression for H
Hs(1)=1;Hs(4)=1; % unit diagonal normalization
assume(Hs,'real');
Vc=diag(sym('Vc', [2 1])); % control regime structural variances
Vp=diag(sym('Vp', [2 1])); % policy regime structural variances
moms=[reshape(Hs*Vc*Hs',4,1);reshape(Hs*Vp*Hs',4,1)]; % compute symbolic expression for moments
moms([3 7])=[]; % eliminate redundant moments
theta=[Hs(2);Hs(3);diag(Vc);diag(Vp)]; % construct parameter vector
D=jacobian(moms,theta); % Take Jacobian
Dt=D;
Dt(1:3,:)=D(1:3,:)*Tc/T; % Scale rows for control regimes moments by share of observations
Dt(4:6,:)=D(4:6,:)*Tp/T; % Scale rows for policy regimes moments by share of observations
grad=matlabFunction(Dt,'Vars',theta); % Convert to numerical function

% testing parameters
alpha=.05; % test level
sig=1-alpha/2; % upper quantile
L=[0 0 -1 0 1 0]; % to conduct wald test of equality of first shock variance across regimes

% create data structures for easy looping
etamat=cat(3,[NF,NY],[RF,NY],[NF,HF],[RF,HF]);
etacmat=cat(3,[NFc,NYc],[RFc,NYc],[NFc,HFc],[RFc,HFc]);
etapmat=cat(3,[NFp,NYp],[RFp,NYp],[NFp,HFp],[RFp,HFp]);

Table5=cell(10,4);
AndrewsTable=NaN(4,4);
%% Estimate each specification

for spec=1:4
    fprintf('Specification %d\n',spec)
    
    etac=etacmat(:,:,spec);
    etap=etapmat(:,:,spec);
    
    etasq1=NaN(Tc,3);
    etasq2=NaN(Tp,3);
    
    etac=etac-mean(etac); % demean so expectation of outer product equals variance (model assumes mean zero in population)
    etap=etap-mean(etap);
    eta=[etac;etap];
    for t=1:T % loop over observations
        v=eta(t,:)'*eta(t,:); % outer product
        if t<=Tc
            etasq1(t,:)=[v(1:2),v(4)];
        else
            etasq2(t-Tc,:)=[v(1:2),v(4)];
        end
    end
    n1=cov(etac,1); % calculate subsample covariances
    n2=cov(etap,1);
    [store,evals]=eig(n2*n1^-1); % compute columns of H as eigenvectors of ratio of var-covar matrices (see e.g. Brunnermeier et al (2019))
    if evals(1)>evals(4)
        store=store(:,[2 1]); % if first variance change is larger than the second, permute columns to match identification assumption
    end
    Hest=[store(:,1)/store(1),store(:,2)/store(4)]; % % normalize for unit diagonal
    S1=cov((Hest^-1*eta(1:Tc,:)')',1); % obtain estimates of structural variances
    S2=cov((Hest^-1*eta(Tc+1:T,:)')',1);
    est=[Hest(2);Hest(3);diag(S1);diag(S2)]; % stack structural parameters into a vector
    [~,~,omega]=momfast(etasq1,etasq2,T,est); % obtain estimate of Omega (moment covariance)
    ghat=grad(est(1),est(2),est(3),est(4),est(5),est(6)); % compute gradient at estimates theta
    V=(ghat'*omega^-1*ghat)^-1; % compute asymptotic variance
    Table5([1 4 7:10],spec)=num2cell(est); % Put estimates in table
    CI1=round([est(1)-tinv(sig,inf)*sqrt(V(1,1))/sqrt(T),est(1)+tinv(sig,inf)*sqrt(V(1,1))/sqrt(T)],4); % CI for H_21
    CI2=round([est(2)-tinv(sig,inf)*sqrt(V(2,2))/sqrt(T),est(2)+tinv(sig,inf)*sqrt(V(2,2))/sqrt(T)],4); % CI for H_12
    Table5(2,spec)={strcat('[',num2str(CI1(1)),',',num2str(CI1(2)),']')};
    Table5(5,spec)={strcat('[',num2str(CI2(1)),',',num2str(CI2(2)),']')};
    if spec>2
        tstat=sqrt(L*est*(L*V/T*L')^-1*est'*L');
        fprintf('Test of equality of first shock variance: t=%.2f\n',tstat)
    end
    % Robust confidence sets
    if spec<3
        range=[-400:.01:400]; % grid for parameter of interest over which to invert test - must be symmetric
    else
        range=[-4:.01:4]; % grid for parameter of interest over which to invert test - must be symmetric
    end
    for pind=1:2
        fprintf('Computing robust confidence set for parameter %d\n',pind)
        [CIbnds, DistCI,stats]=InvertTest(est,etasq1,etasq2,pind,range,alpha);
        Table5(3*pind,spec)={strcat('[',num2str(CIbnds(1)),',',num2str(CIbnds(2)),']')};
        Stats{(spec-1)*2+pind}=stats;
    end
    AndrewsTable(spec,:)=fliplr(min(DistCI(:,1)'>CI2(1),DistCI(:,2)'<CI2(2)));
end
Table4(:,5:8)=AndrewsTable
Table5